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Abstract 

Group theoretic and graphical techniques are used to derive the iV-body wave function for a 
system of identical bosons with general interactions through first-order in a perturbation approach. 
This method is based on the maximal symmetry present at lowest order in a perturbation series in 
inverse spatial dimensions. The symmetric structure at lowest order has a point group isomorphic 
with the Sn group, the symmetric group of N particles, and the resulting perturbation expansion 
of the Hamiltonian is order-by-order invariant under the permutations of the Sn group. This 
invariance under Sn imposes severe symmetry requirements on the tensor blocks needed at each 
order in the perturbation series. We show here that these blocks can be decomposed into a basis 
of binary tensors invariant under Sn- This basis is small (25 terms at first order in the wave 
function), independent of N, and is derived using graphical techniques. This checks the iV 6 scaling 
of these terms at first order by effectively separating the N scaling problem away from the rest of 
the physics. The transformation of each binary tensor to the final normal coordinate basis requires 
the derivation of Clebsch-Gordon coefficients of Sn for arbitrary N. This has been accomplished 
using the group theory of the symmetric group. This achievement results in an analytic solution for 
the wave function, exact through first order, that scales as N° , effectively circumventing intensive 
numerical work. This solution can be systematically improved with further analytic work by going 
to yet higher orders in the perturbation series. 

PACS numbers: 03.65.Ge,02.20,3.75.Hh 
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I. INTRODUCTION 



Obtaining solutions for large systems of interacting particles continues to challenge ex- 
isting approaches and current numerical resources. As the number of particles N increases, 
the Hilbert space that holds an exact solution of the problem scales exponentially in N 
making a direct numerical simulation intractable. [1, 2] For general interparticle interactions, 
this necessitates approximations which in general truncate the Hilbert space of the solution, 
usually by choosing a particular ansatz for the manybody wave function or by truncating 
a perturbation series. Successful methods include mean-field theory and its variants[3], the 
self-consistent multiconfigurational Hartree Fock theory [4], coupled cluster methods which 
are size consistent [5], Monte Carlo methods[6-13], and density functional methods which 
avoid an explicit determination of the iV-body wave function. [14, 15] Perturbation tech- 
niques have typically been based on truncated expansions in the interaction strength and 
thus have been successful in a limited range of the interaction strength. [3] 

In the current paper, we address the challenge of the exponential N scaling for systems 
of N identical bosons by using group theory and by employing graphical techniques. This 
approach is thus inherently analytical in nature although the potential is general. It is based 
on the maximal symmetry present at lowest order in a perturbation series in the inverse 
dimension of space. The lowest-order configuration has a point group isomorphic with the 
symmetric group, Sat [16], with all interparticle distances and angles identical (achievable 
only in higher dimensions). This simple configuration at lowest order yields to an closed 
form analytic solution which can be shown to contain correlation effects. 

The coefficients needed at each order in the perturbation series, can be grouped into 
tensors in the particle label space. These coefficient tensors in the higher-order terms scale 
as (jV 2 )(° rckr+2 ) (N 6 for first order in the wave function) reflecting the N dependence of the 
size of the initial internal coordinate basis. However, due to the expansion about the large- 
dimension, maximally-symmetric structure in which all interparticle distances are equal, the 
size of the internal coordinate basis does not determine the N scaling of the problem. As 
explained in more detail in the body of the paper, the system is order-by-order invariant 
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under the N\ operations of the Sn point group which imposes a highly symmetric structure 
on the tensor blocks needed at each order. We show here that these blocks can be system- 
atically decomposed using a basis of binary tensors invariant under Sjy which are derived 
in this paper using graphical techniques. This basis is finite, independent of N, and in fact 
small (only twenty five terms at first order in the wave function). Since the basis is small, 
there are only a small number of undetermined coefficients in the wave function (twenty five 
at first order) for any iV . These few coefficients are determined by solving the perturbation 
equations. 

Thus by expanding about a large dimension limit which provides the maximally sym- 
metric point group symmetry, we have separated the N scaling problem away from the 
interaction dynamics. Then by allowing the Sn symmetry to do the "heavy lifting", we ar- 
rive at an approach to the iV-body problem which scales as N° . The number of particles N 
enters into the theory as a parameter so that results for a broad range of N can be obtained 
in a single analytic calculation. [17-20] Choosing the perturbation to be in the inverse spatial 
dimension, rather than in the interaction strength, allows this perturbation series to be used 
for both weakly and strongly interacting series as well as the transition between them. 

The final coordinate basis of the problem in which the perturbation equations are solved, 
is a normal mode basis obtained from the zeroth-order perturbation equation which is a 
harmonic equation (a brief four page summary of the method at lowest order, along with 
some lowest-order results, may be found in Ref. 21). The normal modes transform under 
irreducible representations of the Sn group and are obtained using the FG method familiar 
from quantum chemistry. [22] They offer insight into the dominant motions of the iV-body 
system if higher order terms are small. The transformation to normal modes transforms 
the basis of binary invariants into a basis of Clebsch-Gordon coefficients of the symmetric 
group, Sn, which couple the different irreducible representations together to form the scalar 
representation of Sn- This transformation, as reported in this paper, has been achieved 
analytically for arbitrary iV through first order (see Appendices B and C). 

The final result is a solution for the wave function of a system of N identical particles, 
interacting with a general potential and confined in an external potential. Such confined 
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quantum systems are now widespread across many areas of physics and possess from a 
few tens to many millions of particles. Quantum dots [23], Bose-Einstein condensates [24], 
atoms, two-dimensional electronic systems in a Corbino disk geometry [25] and rotating 
superfluid helium systems [26] are all examples of quantum systems confined in an exter- 
nal potential. The interparticle interactions in these systems span an enormous range of 
interaction strengths. The solution involves a large number normal modes, however the 
degeneracy of the corresponding frequencies for a system of iV identical bosons is extensive 
such that there are only five distinct frequencies for a spherical confining potential. Our 
previous analysis shows that these correspond to: breathing and center of mass motions (for 
a harmonic confining potential), single particle excitations (two distinct frequencies) and 
phonon modes. [27] These normal modes will describe the dominant motions of the above- 
mentioned manybody systems if higher-order terms are small. The first-order correction to 
the Hamiltonian has terms cubic in the normal coordinates and their derivatives. This yields 
a first-order, anharmonic correction to the wave function that is equal to the lowest-order 
state times a third-order polynomial of odd powers in the normal modes. 

The approach developed in this paper is generalizable to higher order. This is discussed 
in the Summary and Conclusions. 

The structure of this paper is as follows. In Section II, we set up the formalism for the 
perturbation series in the inverse dimension of space through first order in the wave function. 
Expansions in the inverse dimension of space have been developed for many different systems 
in physics. [28] The approach used in this paper [29, 30] is outlined for a general iV-particle 
problem for L = in this section. In Section III, we discuss the maximal point group 
symmetry of the zeroth order (D — > oo) system and its implications. The presence of 
this extremely high degree of symmetry makes possible the solution of the perturbation 
equations for such a complex, and possibly strongly interacting, system. In Section IV, we 
review the solution for the zeroth-order wave function in preparation for the derivation of the 
first-order, anharmonic correction. In Section V, we develop the mathematical formalism 
needed to solve for the first-order wave function correction and, as an example, apply it 
to deriving the first-order wave function correction for the ground state of a system of N 



5 



identical bosons. In Section VI, we present a summary and our conclusions. 

The formalism developed in this paper, and briefly outlined above, has been tested on 
the example of a three-dimensional fully-interacting confined system of N particles. [31] The 
interaction and confining potential were both chosen to be harmonic. Choosing this sim- 
ple interaction is not necessary or even advantageous for our method, but it did allow an 
independent exact solution of a fully-interacting three-dimensional system for comparison 
with the results of our formalism. This was achieved by comparing the general wave function 
through first-order of Eqs. (33) and (34) with the corresponding first-order wave function de- 
rived by expanding the independent solution to first order. The two, independently derived, 
first-order wave functions are seen to be the same for any N . 

In a second example [32] we show how properties may be derived from the wave function of 
Eqs. (33) and (34), by deriving a closed-form expression for the density profile, an observable 
measurable in the laboratory for systems such as a Bose-Einstein condensate. This general 
result for the density profile is checked by comparing with an independent solution of the 
above mentioned system, N particles interacting harmonically in a harmonic confining po- 
tential. The first-order density profile from the present formalism is indistinguishable from 
the independent solution and shows strong convergence to the exact result to all orders. 

II. THE PERTURBATION FORMALISM IN INVERSE SPATIAL DIMENSIONS 

A. The Dimensional Expansion: Dimensionally-Scaled, Jacobian- Weighted S- 
wave Hamiltonian in Internal Coordinates 

In previous papers we have started the development of the dimensional expansion formal- 
ism for the quantum mechanical problem of N interacting bosons confined by an external 
potential, where N may be less than ten to 10 7 or larger. The zeroth-order wave function, 
density profile and corresponding energy have been previously derived. [17-20, 33] 

In dimensionally-scaled oscillator units the Jacobian-weighted, L = Schrodinger equa- 
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tion reads [19, 34] 
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(2) 



c?7ij »7ifc 

The effective potential V eff is composed of the non- derivative ("centrifugal") portion of the 
kinetic energy U , plus the confining and interparticle potentials, Vconf and V respectively, 
i.e. 

V eif =U + V conf + V, (3) 



where 



N 

^conf v i J 

i=l 



(4) 



at ^N(N-2) + (D-N-l) 2 (^j 



i=l 



(5) 
(6) 



N-l N 

V = Y, E ^(fy). 
i=l j=i+l 

The lim^^oo oc .D 2 to ensure that C/ is finite as D — > oo . Defining <5 to be the inverse 
dimensionality (so that lim^^oo is equivalent to lim^o) an d ((S) to be the ratio of D 2 to 
k(D), 
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then ((5) is finite as 5 — >■ (D — > oo). Defining a position vector, y, consisting of all 
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N(N + l)/2 internal coordinates: 




y 



r 
7 



and r 



, where 7 = 
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\ 7jv-i,jv / 

and taking all of the masses of the particles to be the same, = m , Eq. (1) can be written 

as 

m= (5 2 ((5)f + V eti )<S> = E<S>, (10) 

where 
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V ei M 8, N) = U(y; 5, N) + V cont (y; 5, N) + V(y; 5, N) , 
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(11) 
(12) 

(13) 



dy v = ^- are the derivatives of the elements of the y column vector, and G(y; 5) is the 
N(N + l)/2 x N(N + l)/2 dimensional block-diagonal tensor 
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with I N the N x N dimensional diagonal unit tensor, G 77 (f, 7; 5) the N(N - l)/2 x iV(iV - 
l)/2 dimensional matrix of elements 

[GT"(r, 7; *)](«),<«*) = ^ , (15) 

and 5# is the Kronecker delta (as distinct from the inverse spatial dimension, S). This 
paper uses the implicit summation convention, i.e. repeated indices are summed over (unless 
explicitly noted otherwise), and this convention has been used in Eq. (11). 

B. The Large-Dimension Expansion 

1. The Large- Dimension Limit 

Looking at Eq. (10), we see that as 5 — > (D — > 00) the derivative terms of Eq. (11) drop 
out leaving a static problem where the system localizes at the minima (or more generally 
the extrema) of V e tt(y; 5 — 0, N) at f\ = and 7^ = 7^ . Thus we have the D — > 00 
energy = = [V eff ] gl/2=0 . 

2. Dimensionally-Scaled Internal Displacement Coordinates 

The dimensional expansion is developed by making the following substitutions for all 
radii and angle cosines: 

n = f OD + 5 1 / 2 f' t and 7y=7oo + * 1/2 %-, (16) 
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and denning a displacement vector of the internal displacement coordinates [primed in 
Eqs. (16)] 



V 



and r' = 




where 7' 




(17) 



3. Dimensional Perturbation Series 



After the substitutions of Eqs. (16) and (17) in Eqs. (11)-(15), the Hamiltonian and 
Jacobian-weighted wave function and energy are expanded in powers of 5 l l 2 by deriving 
Taylor expansions in <5 1//2 for G (y(5 x / 2 ); 5) and V ett (y(5 1 / 2 ); 5, N) to obtain 



00 

H = Hoc + 5^ H-i + 5 ^ (6*y Hj 

3=0 

00 

Hr tllij ) =J2(^) J ®j 
3=0 

00 



(18) 



3=0 



The expansion of G (ij(5 l l 2 )\ 5) and V e ± f (y(5 1//2 ); 5, TV)) is most easily achieved using 
the identities of Appendix A. These identities show that to a given order j in <5 1//2 , the 
expansion terms are polynomials of order j of the elements of the y' vector, yl v . If G(y; 5) 
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and V e tt(y;S, N) are functions of 5, rather than 5 1 / 2 , then the equations of Appendix A 
indicate that these polynomials of order j are either all even or all odd in powers of yi v 
depending on whether j is even or odd. 
Thus we have that [35] 

H-i = y' u dy„V eff (y;5 = 0,N)\ y=yoo =0, (19) 
Ho = ~\ dy,d %2 + \ y'„X 2 + ( > , (20) 

Hi = ^G^,^ y' vl dy' V2 d vl z ~ 2 % + ( 3 ^1,^2,^3 ^i^ 2 ^3 + ^ > ( 21 ) 

where the superprescript on the F and G tensors in parentheses denotes the order in 5 1 ^ 2 
that the term enters (harmonic being zeroth order); the subprescripts denote the rank of 
the tensors, and 
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C. The Solution of the Perturbation Equations 

1. The Lowest- Order Wave Function 

The zeroth-order Hamiltonian of Eq. (20) is that of a multi-dimensional harmonic oscil- 
lator. Thus upon transformation to the normal modes of the system, the wave function is a 

product of P = N(N + l)/2 one-dimensional harmonic oscillator wave functions 

p 

$o(y) = II < Mv / ^) , (29) 
11 



where Q v is the frequency of normal mode q' v , and is the oscillator quantum number, 
< n q i v < oo, which counts the number of quanta in this normal mode. 

The transformation to normal modes would appear to be a formidable proposition since if 
N is in the millions, the number of normal modes P is of the order 10 12 or larger. However, 
the D — > oo structure is extremely symmetric (see Section III B) , and it is this maximal 
point group symmetry which allows the normal modes to be derived (see Sec. IV). 

2. Next-Order, Beyond-Harmonic Energy 

The first-order, beyond-harmonic energy correction is 

E 1 = ($ |#i|$o>. (30) 

This is equal to zero which can be seen as follows. The harmonic-order wave function, 
$o ) is a separable product of harmonic oscillator functions. When H 1 is written in terms of 
normal coordinates, it is composed of odd powers of normal coordinates and their derivatives 
(see Eq. (110)). Thus if (f> nv (y/u^q^) is a one-dimensional harmonic oscillator function 

/+oo 
0n„ (V&ql) q„d b ql <Pn„ {V^ql) dq' u = (31) 
■oo 

when a + b = odd, i.e. we have that 

E x = • (32) 

3. Next-Order, Beyond-Harmonic Wave Function 

While the 0(S 3 ^ 2 ) energy correction is zero, the corresponding wave function correction is 
NOT zero. Thus if we write the Jacobian-weighted wave function through first-order beyond 
harmonic as 

$ = (1 + S^A^o + 0(8) , (33) 
then Ai satisfies the commutator equation 

[Ax, H ]Q = H 1 Q . (34) 
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III. THE MAXIMAL POINT GROUP SYMMETRY AND ITS IMPLICATIONS 



A. The Block Structure of the F and G Tensors 



If we generically denote the F and G tensors by Q then we can write the tensors 
in block form as 



(0 2 } Q 



V 



(0) Qrr (0) Qi*7 



\ 



(0) 



For example for the tensor this reads 



(0) 



2 ^vi,V2 



(0)^-y rr (0)^-yr7 
2 2 ^i,(jk) 



\ 



(35) 



(36) 



Likewise the '3 Q tensor may be written in block form as 



(1 3 } Q 



W Qrrr 



W (1) ^ 



(1) Qr-yr 



(1) 



r77 



V 



(1) 



■yrr 



V 



g 77r 



(1) 



7r7 



V 



777 



/ 



/ 



where, for example, the ^ F tensor may be written in block form as 



3 F Vl,U 2 ,V3 



(1) 
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rprrr 
3 r i,j,k 



(!) j^vyr 



( 



y 3^,01),! y 

(1) p7rr 

(!) p77?" 
y 3 r (ij),(M),m y 



(1) j^rry 
3 ^i^(fcl) 

(!) p»*77 
y 3 r i,(jk),(lm) J 

3 (ij),k,(lm) 

(!) p777 
y 3 r (y),(fci),(mn) y 



(37) 



(38) 
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B. The Large- .D Point Group Symmetry, Spj , and the Order-By-Order Invariance 
of the Dimensional Expansion of the Hamiltonian 

The large-.D configuration specified by has the highest degree of symmetry where 
all particles are equidistant from the center of the trap and equiangular from each other (a 
configuration that's only possible in higher dimensions). The point group is isomorphic to Sn 
which in effect interchanges the particles in the large dimension structure. This together with 
the fact that the full D-dimensional Hamiltonian of Eqs. (l)-(8) is invariant under particle 
exchange means the dimensional expansion of Eq. (18) is order-by-order invariant under 
this Sn point group, i.e. the Hj are each invariant under the Sn point group. This greatly 
restricts the F and G tensors of Eqs. (20) and (21). (In three dimensions a corresponding 
iV-particle structure would have a point group of lower symmetry, i.e. one not isomorphic 
to S N despite the fact that all of the particles are identical.) 

C. The Reducibility of the F and G Tensors under Sn 

The maximally symmetric point group Sn , together with the invariance of the full Hamil- 
tonian under particle interchange, requires that the F and G tensors be invariant under the 
interchange of particle labels (the Sn group). In fact the various blocks of the F and G 
tensors are themselves invariant under particle interchange induced by the point group. 
For example, ^ F^ 1 ^ is never transformed into ^^(l/) 7 (fcz)(™™) ' *' e ' ^ e r an d 7 labels are 
preserved since the Sn group does not transform an f[ coordinate into a 7^ coordinate. 

The various r — 7 blocks of the F and G tensors may be decomposed into invariant, and, 
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(39) 



this time irreducible, blocks. Thus for example, may be decomposed into the blocks 

( ?QS v W 
( ?^) v ?<A; ' 

(0 ^S),(.-) v *<^' 

i < j , k < I 

all of which remain disjoint from one another under the group. Significantly, invariance 
under particle interchange requires that tensor elements related by a label permutation 
induced by the point group must be equal. This requirement partitions the set of tensor 
elements for each block into disjoint subsets of identical elements. Consider the elements of 
the ^Q rr block. The element ^QYi belongs to a set of N elements (of the form ^QYi ) 
which are related by a permutation induced by the point group, and therefore must have 
equal values. Likewise, the element ^Q"^ belongs to a set of N — 1 elements related by a 
permutation induced by the point group and sharing a common value. Proceeding in this 
fashion, we observe that the blocks of the harmonic-order tensors are partitioned into the 
following set of identical elements which remain disjoint under particle interchange: 

V i and k (40) 





- {0) Q 

— 2 V, 


(O)^rr 
2 <Aj 


= «$Q, 


(0)^r7 
2 **i,(ik) 


= 


(0) n r 7 
2 Vi,(,-fc) 


= ( s>q: 


(0) ^77 


- (0) O' 

— 2 V 


( )/O77 
2 ^(ij),(ik) 


= °$Q 


(0)/O77 
2 ^(ij),(kl) 


= °$Q 



ft V ^jandfc^ (41) 

Sim) = { °2QZpn) V » < ^ > Km, a ndn>p (42) 

IJmn) V j <k and I ^m<n (43) 

$).(«) V i<jandfc<Z (44) 

JiL),(in) V i < j ^ k> i and I < m ^ n > I (45) 

5n),H V i*j*k*l andm^ntptq (46) 
i < j , k < I m < n , p < g 
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The harmonic-order block matrices contain the sets of elements in Eqs. (40)- (46) arranged 
in an intricate pattern. In Ref. 17 it was shown that this arrangement could be expressed 
in terms of binary matrices formed from matrices commonly used in graph theory. In this 
paper, we introduce a generalized structural decomposition for higher orders by introducing 
binary tensors, derived using graphs, i.e. an extension of binary matrices. 

D. Introducing Graphs 

Definition 1 A graph Q = (V,E) is a set of vertices V and edges E. Each edge has one 
or two associated vertices, which are called its endpoints. 

For example, J is a graph Q with three vertices (or "dots") and two edges (or lines). 
We allow our graphs to include loops and multiple edges [36]. A graph contains information 
regarding the connectivity of edges and vertices only: the orientation of edges and vertices 
is insignificant. Two graphs with the same number of vertices and edges that are connected 
the same way are called isomorphic. 

We introduce a mapping which associates each tensor element with a graph as follows: 

1. draw a labeled vertex (• i) for each distinct index in the set of indices of the element 

2. draw an edge ( l_J ) for each double index {if) 

3. draw a "loop" edge ( Q) i ) for each distinct single index i 

For example, the graph corresponding to the tensor element ^Q^ij) under this mapping is 

CT 3. 

This mapping will be used to represent both the value of sets of identical tensor elements 
and their tensor structure. 

The partition of each set of tensor blocks is composed of a certain number of sets of 
identical elements. Under this mapping, the image of each set of identical elements is a 
set of isomorphic graphs. We represent each set of isomorphic vertex-labeled graphs (and 
therefore the corresponding set of identical tensor elements) by an unlabeled graph. 
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We denote the set of unlabeled graphs for each block as &x 1 x 2 ...x R , where R is the rank 
of the tensor block (and therefore the number of edges in each graph in the set). 

G rr = {§,0 0} 

G 7 r = {O-, O —} (47) 

g 77 = {o, j,~} 

G r = {0} 
G 7 = {-} 

G rrr = {Q^,gO,0 0} (48) 
G 7 , r = {^,00,0-0,00,00} 

c 77r = {oo>oo»^.l-o,LO.^»no} 

These graphs provide a convenient notation to label the sets of tensor elements. We represent 
the value of a set of tensor elements by the corresponding unlabeled graph for that set. For 
example, the value of the tensor elements ^QT^j) * s represented as ^QiO-^)- We list 
the graphs and corresponding elements in Tables I, II and III. Table I lists the graphs for 
matrices at harmonic order, along with the tensor elements represented. 

At the first anharmonic order we have the graphs and the corresponding rank-one and 
rank-three Q-tensors in Tables II and III respectively. 

In the next subsection, we represent the tensor structure of each set of identical elements 
corresponding to a graph Q with the binary tensor B{Q). 

E. Decomposition of Q in the Basis of Binary Invariants 

Each tensor block is composed of disjoint sets of identical elements, each labeled by a 
graph. Each block has a finite number of such sets (and corresponding graphs). We exploit 
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these symmetry properties by introducing a binary tensor for each block and for each graph 
that embodies the structural arrangement of the corresponding elements. 

Definition 2 A binary invariant [B block {Q)] vltV2} „, corresponding to a graph Q is a tensor in 
which elements related to Q by the above mapping are unity and all other elements are zero. 

A binary invariant B bloch (Q) of a tensor ^ der ^Q has the following properties: 

• The rank of the tensor B block (Q) equals the number of edges of the graph Q. 

• B block {g) and B block (Q') are linearly independent, since binary invariants of different 
graphs share no common elements. 

• The set of binary invariants for the graphs of a block span the unit matrix (that is the 
sum of a set of binary invariants yields the unit matrix). 

• B block (Q) is invariant under particle interchange. 

We introduce a method for deriving expressions for these binary invariants and list them 
through first-anharmonic order in Ref. 37. 

Therefore, we may resolve the Q tensors at any order as a finite linear combination of 
binary invariants: 

'r'QIuL,^ = E Q UOCk W [B block (G)] uu ^ R (49) 

Q&Gbloak 

where <&biock represents the set of graphs present in the order-O, rank-i? tensor block 
^Q block , and the binary invariant B block (Q) has the same dimensions as the original Q 
tensor block. The scalar quantity Q block (Q) is the expansion coefficient. 

The resolution of symmetric tensor blocks in the basis of binary invariants in Eq. (49) 
represents a generalization of a technique used at harmonic order in Ref. 17 to arbitrary 
order. This equation also separates the specific interaction dynamics present in Q block (Q) 
from the point group symmetry embodied in B block (Q) . 
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IV. LOWEST ORDER: NORMAL MODES AND THE HARMONIC WAVE 
FUNCTION 



According to Sections II B 1 and II C 1, the Jacobian- weighted wave function of spherically 
confined quantum systems in large dimensions becomes harmonic, with the system oscillating 
about the D — > oo configuration where every particle is equidistant and equiangular from 
every other particle. Notwithstanding the relatively simple form of the large-dimension, 
zeroth-order wave function, it includes beyond-mean-field effects. 

^From Eq. (20) one sees that the Schrodinger equation for the zeroth-order Jacobian- 
weighted wave function has the form 

[~\ d,,dy, 2 + \ ^F VuV2 y' v X 2 + ( > ) $ (y') = E $ (y') (50) 

where y' is the dimensionally-scaled displacement coordinate vector. The matrices (i.e. rank- 
two tensors) and in Eq. (50) are both symmetric constant matrices, while 
is a scalar quantity. 



A. The F G Method 

The Wilson FG method[17, 22] is used to solve equation (50) for the normal mode coor- 
dinates and frequencies and allows us to exploit the point group symmetry of the problem. 
The 6 th normal mode coordinate, q' b , may be written: 

ql = V b ,„y' l/ . (51) 

The coefficient matrix V satisfies the eigenvalue equation 



(o) F (o) G 



[V T U, b = Kb) [V T ] Uub (52) 



and orthonormal condition 

Sab ■ (53) 
Eigenvalues, \(b) , are related to the frequencies, oj{b) , by 

A(6) = co(b) 2 . (54) 
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In principle, equation (52) is still a formidable equation to solve unless iV is quite small since 
there are P = N(N + l)/2 normal coordinates and up to P distinct frequencies. 

B. Group Theory: The Sn Symmetry at Harmonic Order 

As defined in Eq. (16), the large-dimension structure is completely symmetric, resulting 
in an Sn point group. Also, as we have noted above, the full S'-wave Hamiltonian is invariant 
under particle interchange, an Sn symmetry under which the system is invariant. These two 
facts mean that Eq. (50) is also invariant under the group Sn- This allows us to use the 
group theory associated with Sn symmetry which brings about a remarkable reduction from 
P possible distinct frequencies to five actual distinct frequencies. The Sn symmetry also 
greatly simplifies the determination of the normal coordinates. The vehicle for this is the 
use of symmetry coordinates (see below). [22] 

The fact that , ^G are invariant matrices under S N , implies that the ^F ^G 
matrix of Eq. (52) is also invariant under Sn ■ This means that the ( ° } F ( ° } G matrix may 
be expanded in terms of the second rank binary invariants. The ^F , ^G , and ^F ^G 
matrices, which we have generically denote by > are P x P matrices which may be 
written: 



The upper left block of Eq. (55) is an N x N matrix with elements denoted by indices i,j . 
Defining the number of 7^ coordinates to be M = iV(iV — l)/2 , the upper right block is an 
N x M matrix with elements denoted by indices i, (jk) . The lower left block is an M x iV 
matrix with elements denoted by indices (ij), k . The lower right block is an M x M matrix 
with elements denoted by indices (ij) , (hi) . Using Eq. (49) we can write 




(55) 




Ve[B(0-)] T + /[ J B( ~)] T gB( )+hBU) + LB(~)J 




(56) 
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where a, b c, d, e, f , g, h, and i are given by Eq. (42) of Ref. 17. The S jv point 
group symmetry means that the ^F ^G matrix, instead of having a possible P x P = 
N(N + 1)/2 x N(N + 1)/2 independent elements, in fact requires that only nine independent 
quantities have to be calculated! 

The fact that the ^F ^G matrix of Eq. (52) is an invariant matrix under Sn further 
implies that the eigenvectors V of the matrix and the normal modes transform 

under irreducible representations of SV Using the theory of group characters [19] one sees 
that the are reducible to one one-dimensional irreducible representation labeled by [N] , 
and one (N — l)-dimensional irreducible representation labeled by [N — 1, 1] . Here, the 
partition, [oi, a 2 , a 3 , ... , a m ] where a\ + a 2 + a 3 + ■ ■ ■ + a m = N , denotes the irreducible 
representation of .[16, 38] The 7^ are reducible to one one- dimensional irreducible rep- 
resentation labeled by [TV] , one (N — l)-dimensional irreducible representation labeled by 
[N— 1, 1] , and one N(N— 3)/2-dimensional irreducible representation labeled by [N—2, 2] . 
Since the normal modes transform under irreducible representations of and are composed 
of linear combinations of the elements of the internal coordinate displacement vectors f' and 
7' , there will be two 1-dimensional irreducible representations labeled by the partition [N] , 
two (N — l)-dimensional irreducible representations labeled by the partition [N — 1, 1] 
and one entirely angular iV(iV — 3)/2-dimensional irreducible representation labeled by the 
partition [N — 2, 2]. All of the normal modes transforming together under the same irre- 
ducible representation will have the same frequency, and so rather than N(N + l)/2 distinct 
frequencies there are in fact only five distinct frequencies! 

C. The Program to Determine the Normal Modes 

We determine the normal coordinates and distinct frequencies in a two-step process: [19] 

a) We start with an important observation: both f ' and 7' transform under or- 
thogonal and, as we have seen above reducible representations of the symmetric 
group, since both f' T . f' and 7 /T .7' are invariant under any permutation of par- 
ticles. Thus y' also transforms under an orthogonal reducible representation of 

21 



the Sn group. 

Then we determine two sets of linear combinations of the elements of coordinate 
vector f' which transform under particular orthogonal [N] and [N — 1, 1] irre- 
ducible representations of Sn ■ These are the symmetry coordinates for the r 
block of the problem [22]. Using these two sets of coordinates we then determine 
two sets of linear combinations of the elements of coordinate vector ~f which 
transform under exactly the same orthogonal [N] and [N — 1, 1] irreducible rep- 
resentations of S N as the coordinate sets in the r block. Then another set of linear 
combinations of the elements of coordinate vector j' , which transforms under 
a particular orthogonal [N — 2, 2] irreducible representation of S N , is derived. 
These are then the symmetry coordinates for the 7 block of the problem[22]. 
Furthermore, we choose one of the symmetry coordinates to have the simplest 
functional form possible under the requirement that it transforms irreducibly 
under Sn ■ The succeeding symmetry coordinate is then chosen to have the next 
simplest functional form possible under the requirement that it transforms irre- 
ducibly under Sn , and so on. In this way the complexity of the functional form 
of the symmetry coordinates is kept to a minimum and only builds up slowly as 
more symmetry coordinates of a given species are considered. (More on Step a) 
in Section IV D.) 

b) The matrix originally expressed in the r', 7' basis is now expressed 

in the symmetry coordinate basis resulting in a remarkable simplification. The 
N(N + l)/2 x N(N + l)/2 eigenvalue equation of Eq. (52) is reduced to one 2x2 
eigenvalue equation for the [N] sector, N — 1 identical 2x2 eigenvalue equations 
for the [N — 1, 1] sector and N(N — 3)/2 identical lxl eigenvalue equations for 
the [N - 2, 2] sector. For the [N] and [N — 1, 1] sectors, the 2 x 2 structure of 
the eigenvalue equations allows for mixing of the symmetry coordinates in the 
normal coordinates in the r and 7 blocks. The lxl structure of the eigenvalue 
equations in the [N — 2, 2] sector reflects the fact that there are no [N — 2, 2] 
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symmetry coordinates in the r block to couple with, i.e. the [N — 2, 2] normal 
modes are entirely angular. (More on step b) in Section IV E.) 



D. Transformation to Symmetry Coordinates and the W Matrices. 



Symmetry coordinate vectors S° and of the r and 7 blocks respectively, where a and 
(5 are the partitions labeling the irreducible representations of Sn under which the symmetry 
coordinates transform, are said to belong to the same species when a = (5. In what follows 
we will find it useful to define a full symmetry coordinate vector as follows: 

/ Sl N] \ 



s = 




( 



S [N] 
g[N-l, 1] 
S [N-2, 2] 



\ 



= Wy', 



(57) 



where 



S [N] 



N] 



Si: 



S [N-1, 1] 



and W is the orthogonal matrix 



/ clN-l, 1] 
q [N-l, 1] 



( wl N] 





and 2 1 = - 2 " 2 1 , 



(58) 



W 







w. 



[N] 



Wr 



[N-l, 1] 









[N-l, 1] 



V 







[N-2, 2] 



(59) 



In S, symmetry coordinates of the same species are grouped together. The Wx> , where 
a = [N] , [N - 1, 1] or [N - 2, 2] , and X' is r or 7 (only 7 for the [N - 2, 2] sector), are 
given by 

[Wl% = ^=[lrV and [Wf]] fe) = ,/— - ^— [l 7 ] (y) (60) 
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for the [N] species, where 



[l r ]i = l Vl<i<iV and [1 7 ] W) = 1 Vl<i,j<N. (61) 
Regarding the [N — 1, 1] species 

[Wf" 1 ' 1] k= ^ fc^- j W = /T i- r r (Q»-fc + i-^i + i,*) . (62) 
where 1 < i < N - 1 and 1 < k < N , 



(63) 



m=l 

= when j — i > 1 , 

and 

[wf - 1 ' \ m = 7 ^ ([w^> % k [i r } l + [w^> ^uirU) 

= ~TT- 1 WAr =HT ( (©i-fc+1 [lr]i + ©i-I+l [lr]fc) ~ « (<*i+l, k M + <&i+l, i Wfe) J ,(64) 

0(i + 1)(JV - 2) V / 

where [wj^ -1 ' ^} ik and [l r ]j are given by Eqs. (62) and (61) respectively, and 1 < k < I < N 
and 1 < % < N - 1 . Finally, for the [TV - 2, 2] sector 

[w\ N - 2 ' % ]Umn) = .... = = — ( (e,_ m+1 - ^ +1 , m )(e,_„ - (j - 3)^- B )+ 

0(1 + 1)0-3)0-2) v . 

+ (©i-n+l - t5i+l, n )(Qj-m ~ (j - 3)6 jm )\ , 

(65) 

where l<i<j-2,A<j<N and 1 < m, n < N . 



E. Solution of the W- Transformed Eigenequations for the Normal Coordinates 

The symmetry coordinate transformed ^F W ^G w = W ^ ^F ^G j W T matrix re- 
mains invariant under Sn and so may be expanded in terms of the Clebsch- Gordon coef- 
ficients of Sn coupling two irreducible representations a and (5 together to form a scalar 
(i.e. invariant) [N] representation, i.e. these Clebsch-Gordon coefficients are invariants in 
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the symmetry coordinate basis. According to group theory [39], the only non-zero bilinear 
invariants which couple two vectors transforming under orthogonal (more generally unitary) 
irreducible representations, couple together equivalent irreducible representations. Further- 
more, there is only one non-zero bilinear invariant which couples two vectors transforming 
under equivalent, orthogonal irreducible representations. As the irreducible representation 
is orthogonal, this invariant is proportional to the Kronecker delta (unit matrix), J a .[40] 
Thus we can write 



(o) F (o) r 

2 W 2 ^ W 
( (P)„FG 



2 a [N][N] ® I[N] 



V 



o 
o 



o 

(0)_FG 

2 °[N-1, l][N-l, 1] 





I 



[N-l, 1] 








\ 



(0) FG 



If 



(66) 



2 U [N~2, 2][N-2, 2] ^ A [N-2, 2] J 

We then write Eq. (52) in the symmetry coordinate basis of Eqs. (57)-(65) where it reduces 



to three eigenvalue equations of the form 



(o) fg" 

2 u aa 



[c a ] YXl = \(a,Y) \c a } YX2 



The reduced (lJ 'G block matrices 



(0), 

2 



(0) a FG 



and 



(0) „FG 



2 °\N- 1,1] [AT- 1,1] 



(67) 



are 2x2- 



- X2,X\ 



(0)_FG 

2 u [7V-2,2][7V-2,2] 



is a 



X2,X\ 

(0) fg" 



2 u [Af][Af] 

dimensional matrices (Y — ± and X x , X 2 = r or 7) , while 

one-dimensional matrix (Y = 7 and Xi — X 2 — 7). The elements of the 
matrices are known analytic functions specific to the individual systems being considered. 
Thus there are five solutions to Eq. (67) which we denote as = {A([iV],±), [c^]±} , 
1± = {A([A^ - 1, 1], ±), [c^" 1 ' 1] ]±} and 2 = {A([iV - 2, 2]), c^- 2 ' 2 ]} . The two element 
[c a ]± vectors for the a = [N] and [N — 1, 1] sectors determine the amount of angular-radial 
mixing between the symmetry coordinates in a normal coordinate of a particular a . Thus 



[q'] b = [c«] ±r [S?fc + [c«] ±7 [S«] { 



(68) 



The normal-coordinate label, 6, has been replaced by the labels a, £ and ± on the right 
hand side of Eq. (68). 
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^From Eqs. (53) the c a also satisfy the orthonormal condition 



[c a ]yiXi [c a ]y 2 x 2 



(0) G 
2 °aa 



X2,Xl 



2 a aa 



>Y U Y 2 , 



matrices are the reduced matrices in the symmetry 



(69) 



where the diagonal 
coordinate basis. 

For the [N — 2, 2] sector, the symmetry coordinates are also the normal coordinates up 
to a normalization constant, c^ -2, 2 ' (see Eq. (69) above). 



F. The Wave Function. 

Thus the wave function of Eq. (29) is the product of P = N(N + 1)/2 harmonic oscillator 
wave functions [19] 

$o(y')= II ( 7 °) 

At=0±,l±,2 

where \i labels the manifold of normal modes with the same frequency u) M and 

M%) = U^AV^<i'J » ( 71 ) 
i=i 

0n M iV^j.^i) 1S a one-dimensional harmonic-oscillator wave function of frequency uj^ and 
is the oscillator quantum number, < n m < oo, which counts the number of quanta 
in each normal mode. The quantity is the number of normal modes in the degenerate 
manifold of frequency frequency , and has values = 1 , N — 1 or N(N — 3)/2 for 
pi — O 1 * 1 , l 1 * 1 or 2 respectively. 



V. NEXT-ORDER: THE FIRST ANHARMONIC WAVE FUNCTION 

As we have seen in Section II C 2, the next-order, beyond-harmonic energy correction is 
zero. This though, is not the case for the next-order, beyond-harmonic wave function. 
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A. Next-Order, Beyond-Harmonic Wave Function 

Expanding the wave function through first-order beyond harmonic according to Eq. (33), 
the wave function correction Ai, satisfies the commutator equation Eq. (34). As we have 
noted, the lowest-order wave function $o(y') (see Eqs. (70) and (71)) and Hamiltonian H 
are separable in the normal modes q' . Thus to most easily calculate the effect that the 
next-order Hamiltonian H x (see Eq. (21)) has on the wave function (and energy) we need 
to express Hi in terms of the normal coordinates. 

As in the discussion of the harmonic-order solution, we achieve this in a two step process, 
first transforming to symmetry coordinates (an orthogonal transformation) and then to 
normal coordinates. 



B. Transformation of Hi to Symmetry Coordinates 

/From Eqs. (57) and (59) 

y' ui = W^ U2 S V2 (72) 

and 

dvW = ^ > ( 73 ) 

where we remind the reader that we are using the repeated index summation convention 
and 1 < ui, v 2 < P = N(N + l)/2 , then from Eq. (21) 



1 rWc...! c c c , K 1 ) 



where 



[?G w l = W v , v [?G] V (75) 

[?F W ] V = W vv ^ ] F] V (76) 

[3 ^w\v\,Vl,VZ — Wi/l,JJl ^V2,V2 W U3 ,r]3 [3 ^}vi,V2,V3 (77) 

[3 Fw]vi,V2,V'i — ^V2,m Wi/ 3jJ]3 [3 F}rn,r]2,m ' (78) 
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and 1 < 771,772,^3 < P ■ I n Sections III B and III C we noted the fact that the F and G 
tensors are invariant under Sn, leading to the expansion of these tensors in terms of binary 
invariants B{Q) specified by the graph Q (see Eq. (49)). Thus the transformation of the F 
and G tensors results from the transformation properties of the binary invariants. 
The third-rank binary invariants transform as 

W„ um W U2 , m W vzm [B{G)\ (79) 

while the rank-one binary invariants transform as 

[B w (Q)]u = w v , v [B(g)] v (80) 

Now if T is an element of the S N group such that 

y' T = Ty', (81) 

then the symmetry coordinate vector S transforms as 

S T = T W S (82) 

where 

T W = WTW T . (83) 

It follows that the symmetry coordinate transformed binary invariant B W {G) is also invari- 
ant, coupling three symmetry coordinates, or one symmetry coordinate, together to give a 
scalar term under Sn, i-e. one that transforms under the [N] representation. Since the S 
vector is composed of symmetry coordinates transforming under irreducible representations 
of Sn (see Eqs. (57) and (58)), B W {Q) is thus composed of Clebsch-Gordon coefficients of Sn 
which couple the irreducible representations together to form a scalar [N] irrep. Regarding 
the rank 1 terms, where there is only one symmetry coordinate involved, this has to be Sy^ 
or since the other symmetry coordinates are not scalar under Sn- Thus we can write 
in the block form of Eq. (57) 

B W (0)= ([iV(O)] r ,0,0,0,0) , (84) 

28 



and 



B 



w { 



0, 



ay 



,0,0,0 



where from Appendix C, Table V 



iV(O) 



N. 



5V 



N(N-l) 



so that 



and 



b w (o)s= SV(o) 



B W (~)S = 



Si* 



s[ N K 



(85) 

(86) 
(87) 

(88) 
(89) 



Now consider the rank 3 terms. There are only so many ways we can couple three 
irreducible representations drawn from [N] , [N— 1, 1] , and [N — 2, 2] irreps. to form a scalar 
[N] irrep. . From the Clebsch-Gordon series for S N by Murnaghan, Gamba and others [41], 
these are 









[7VH 


5 [AT] $ 


5[JV] 


= [N} + 






5 [-/V — 


1,1] e 


§ [AT — 


1,1] 


= [N] + 




[N}$ 


§ [AT — 


2,2]$ 


§ [AT — 


2,2] 


W] + 


[N- 


1,1] 


5 [AT — 


1,1] 


5 [N — 


1,1] 


= W] + 


[N- 


1,1] 


§ [N- 


1,1] 


3 [AT — 


2,2] 


= [N} + 


[N- 


1,1] 


5 [AT — 


2,2]$ 


3 [Af — 


2,2] 


= [N] + 



(90) 



and two linearly independent couplings 

[N - 2, 2] <g> [N - 2, 2] <g> [A/ - - 2, 2] = 2 [AT] + 



(91) 



In what follows we denote these two different couplings of three [N — 2, 2] irreps. together 
to form an [N] irrep. by the roman numerals for the numbers 1 and 2 , i.e. J and 11 
respectively. Thus we can write a third rank Bw(G) in block form as 
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B W {Q) 



( 



B%°(G) = o 



I 



\ 



K B™(g) = o / 



( 



BT{<3) = o 



\ 



B° W 12 (G) = 



\ 



( 



\ 



B^\G) 

\ B%°{G) = ) 
( \ 



B%°(G) = 



^ B 2 W 20 (G) ) 



\ 



BW\G) 
B%\G) 
K B^(G) / 

( x 

B$\G) = 



B° W 2 \G) 



B$\G) = 



B^\G) 



\ 



B^\G) 



( \ 

B1§\G) 



J 



B 2 W 12 (G) 



B 2 W 22 (G) 



(92) 



/ 



where each block B^ a2a3 (G) is related to the appropriate Clebsch-Gordan coefficient by a 
scalar multiplier /3 aia2a3 ^' R -(Q^ 1 



aia2Ce 3 ,TZ 



(93) 



n 



In Eq. (93), C aia2a3 ' n is the appropriate Clebsch- Gordon coefficient of Sn for the irreps. . 
As indicated above there are two linearly independent Clebsch- Gordon coefficients coupling 
three [N — 2, 2] irreps. together to form a scalar [N] irrep. Thus 1Z has two values, running 
from / to II when a X 0L 2 a z = 222 . All other couplings of the [N] , [TV — 1, 1] and [N - 2, 2] 
irreps. have only one Clebsch-Gordon coefficient, and so 1Z has only one value, I . 
One further notes that 

C aia * a3 > n = ( ai a 2 ) C aia2as ' n (94) 
where (aia 2 ) is the operator which interchanges a.\ with a 2 . Likewise for all of the other 
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permutations of a>i, a 2) and a 3 . Thus we have 



(95) 



and likewise for the other permutations of a\,a 2 , and a 3 , i.e. f3 a ^ a 'i< n {g) is symmetric in 
cti, CK2, and «3. 

The necessary Clebsch- Gordon coefficients are derived in Appendix B, while the multipli- 
ers, f3 OLlOL2az ' 1l {Q) , of the Clebsch-Gordon coefficients in the expansion of the binary invariants 
in the normal coordinate basis are derived as analytic functions of N in Appendix C. These 
calculated values, together with the values of the first anharmonic G and F elements for the 
particular system of interest, result in a first-anharmonic order Hamiltonian in the symmetry 
coordinate basis. 

Thus pulling everything together from Eqs. (49), (74)- (80), (84), (85), and (93), we find 

[ 3 ^Wlviww ~ 



G 



n 



, (96) 



[Wp l 

[ 3 r W\vt,V2,vz 



Xi,X2,Xs 



^1011012013,12. 



, (97) 



[?G W ] U = E {l iG x {G) [Bfr{Q)]t 



[?F W 1 = E { ^F X {Q) [B$,t9)\t = 
g&Gx 



(1)_G 
1 °ot 



1 °ot 



c?, 



(98) 
(99) 
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where 



'(1)_G 

3 OLlOL20LZ,Tl 



(1)_F 

3 ^aict2C»3,7?- 



Xl,X2,X3 



Xl,X2,X3 



(1) (7 G 
1 °a 



1 °a 



GeG Xl x 2 x 3 



Xl,X2,X3 



, (100) 



X\,X2,Xz 



101) 



£ @) [IV (<?)] 



(102) 



We also have 



Cf = 1 when ct = , and = when a = 1 , or 2 , 



(103) 



since, from Eqs. (84) and (85), Bfo(G) = when a = 1 or 2 . In Eqs. (96)-(101) above, 
repeated ol\ , oli , 0:3 , Xi , X2 , X3 , and X labels are not summed over. We remind the reader 
that X 1 , X 2 , X 3 , and X label r or 7 ; a t , a 2 , «3 , and a label the irreducible representation 
of Sjv ([N] = , [X - 1, 1] = 1 , or [N - 2, 2] = 2). The labels v 1 ,u 2 ,v 3 , and z/ cover 
both eti , ol 2 , (X3 , a , and Xx , X 2 , X 3 , X as well as the irrep. element labels i , j , , etc. . 
A summary of our use of different label sets can be found in Ref. 37. 



C. Transformation To Normal Coordinates 

The transformation to normal coordinates is affected by a transformation which only 
mixes the r and 7 blocks i.e. 

q a = C a S a (104) 

where S a is given by Eq. (58). The matrices C° and C 1 are 2x2 matrices, while C 2 is a 
scalar quantity. 
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For the and 1 sectors 



Thus 



and 



Defining 



/ 



C a = 



\ 



C°cos0£ 
C?cos0° 



\ 



C% sin 9 
Ct sin 6 



J 



S a = {C a )- 1 q a 



d 



d 



dS a dq a 



(105) 



(106) 



(107) 



1 'a 
1 'a 



J Y 



= c 



1 OL 



"(1) T G 
"(l) r F 



= [C 



-llT 



J x 

(1) F 
1 °a 



\l) a G 

3 ai«2a3,1? 



- P 1 ]Yi,Xi^V2,X 2 Cy3,X3 



(108) 



3 aia^a^,^ 



^1 1^2,^3 



and defining Gy and Fy to be the G and F tensors in the normal coordinate basis, 



i "^v 
1 b v 



1 r a 



1 'a 



(1) 



3 



3 



f 1,^2,^3 



n 



y 

(1) T G 

3 aiQ2Q3,"R 



(109) 



3 aia2<^3 ? ^- 



n,y 2 ,y 3 



n,y 2 ,y 3 



]?1,6,«3 ' 



r / ^aia203>^-l 



where repeated a± , a 2 , and a 3 labels are not summed over, we obtain H 1 in normal coor- 
dinates 



1 



-^Gvld,,+[^F v ]„ql, 



(i) 



(110) 



/From Eqs. (33) and (34), 
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$ = (l + ( J 1 /2 Al )$ + O ( ( j) > (33) 

[Ax, # ]$o = #i$o, (34) 
and (110), we can derive the first-order, beyond-harmonic wave function. 

D. An Example: The First-Order, Beyond-Harmonic Wave Function for the N- 
boson Ground State 



Now 



and 



dq 



#<l>o(u)' q') = u(uq' 2 - l)(j>o(u' q ) . 



Ill) 



(112) 



Therefore from Eq. (110) we find that the action of Hi on the lowest-order, ground-state 
wave function g $ becomes equivalent to the action of a 3rd-order polynomial (-^i) eff : 



Hi g $ = {Hi) cS 3 $ 



(113) 



where 



(ft) 



off 



1 

2 



3 U V 



(l) r 
3 



^2^3 + T7T 

Vl,V2,V3 OI 



(1) F 

3 



^2 + 2 



(1) 



v 



VI, V2, v-i 
V V1 + 



i 



(114) 



We use the index /x to refer to both the partition a and the block Y: /i e {0+, 0+, 1+, 1—, 2}. 
This is a convenient index to refer to the 5 blocks of the normal mode vector q'. In this 
notation, the elements of the block q'^ are denoted [q'^\^, with the range of £ implied by /i. 
Using this notation, we define the tensor blocks ^t^ 1 n and ^V^ 1 : 



(l) T Hi 

3 m,H2,H3,T^ 

1 r o± 



1(1) G J_ (1) F 

2 3 T //l, J U 2 ,/i3,7?. a; M2 a; M3 + J) 3 T Ml,/i2,^3,^ 



(1) t G r j -I- (1) t G m , 4- (1) <r F 
3 T 0±,m,M W M + 2 1 ^i^Oi + i T 0± 



(115) 
(116) 
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where d ± = 1 , d 1± = N - 1, and d 2 = N(N - 3)/2. 
/From Eqs. (34) and (114) 



Ai= E ( 117 ) 



where 



Ml>M2,M3 

71 



_ (l) T Hi 

(1) A __ 3 /ji,/j 2 ,A t 3,^ 



J) A _ J Ml»^2,M3,/< /110\ 

(!) T A -IfjiU+Vd f (1) r A + (1) r A + {1) r A M (119) 

1 r 0± — ^ 1 r 0± + 2.^1 V\ 3 r 0±^ + 3 V0±/i + 3 V^tOi y ^ l iiy J 

Notice that the next-order wave function is equal to 9 $o times a polynomial of odd powers 
of order three. In Eqs. (115)-(119) repeated indices are not summed over unless explicitly 
indicated. 



VI. SUMMARY AND CONCLUSIONS 



The iV-body problem is ubiquitous in physics as well as chemistry, biology and related 
fields. Due to the exponential N scaling of these problems, their solutions require approxi- 
mations. Our work suggests that it may be advantageous to compartmentalize the N scaling 
problem away from the rest of the physics. We have shown that this may be done for a 
system of N identical bosons, using group theory and graphical methods. The N scaling is 
checked due to the structure imposed at each order of the perturbation series from the invari- 
ance under the large-dimension point group. Having identical particles is a necessary, but 
not sufficient condition to take full advantage of the group theory of the symmetric group. 
The interparticle distances and angles must be identical, i.e. a point group achievable only 
in higher dimensions. [17, 33] This point group in high spatial dimensions is isomorphic with 
the symmetric group, Sn- 
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As we have seen above in the body of this paper, higher-order terms involve a polynomial 
in coordinates quantifying the deviation of the system from this lowest-order configuration. 
The coefficients weighting individual monomial terms in these polynomials may be gathered 
together to form tensors in the particle label space, where the particle label i runs over 
1 < i < N . Since these tensors are invariant under the N\ operations of the Sn group, they 
yield to an exact decomposition in terms of a small basis of binary tensors, invariant under 
Sn. The basis of binary invariants is finite, independent of N, and has been derived using 
graphical techniques. Thus the problem and its solution involve a small, iV-independent, 
number of coefficients for this basis. We thus now have a problem which scales as iV° where 
N enters as a parameter. 

The result of the analysis outlined in this paper is the wave function for a general, inter- 
acting system, exact and analytic through first order. This first-order correction introduces 
anhamonicities and couples the breathing, single-particle and phonon mode coordinates to- 
gether, which for the ground state takes the form of a cubic polynomial of odd powers of the 
normal mode coordinates multiplying the lowest order solution (see Eqs. (33) and (117)). 

In principle any observable quantity of a confined quantum system can be derived from 
this correlated wave function. The density profile of a macroscopic confined quantum system 
offers an experimentally accessible window into beyond-mean-field effects in a quantum 
object. These effects include a greater spatial extent resulting from hard collisions within 
the system as well as fermionization and crystallization in quasi-one and two dimensional 
systems [42]. Results for the density profile through next-order for large N and/or large 
interparticle interaction strength may be derived. [32] One notes that the derived Clebsch- 
Gordon coefficients couple up to three normal modes together to form a scalar [N] function 
and this will allow excited-state wave functions of up to three quanta to be derived for 
systems which have completely spatially symmetric wave functions, i.e. as a spinless bosonic 
system (e.g. a BEC). The fact that excited states may be derived means that thermal 
properties at finite temperatures for cold confined quantum systems can be calculated in an 
appropriate thermal sum over states. 

This method generalizes to higher order. Successive terms in the perturbation series will 
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require a finite number of binary invariants as a basis for the tensor blocks at each order. 
The number of binary invariants at a particular order is independent of N, thus checking 
the N scaling. However, there is still significant analytical work to be done in deriving 
the Clebsch- Gordon coefficients at each order. This work is only possible through the use 
of the group theory of the symmetric group which greatly constrains the structure of the 
problem and restricts the number of coefficients. Although the work needed will increase 
at each order, it remains analytic and will be a systematic extension of the work outlined 
in this paper for the first-order correction to the wave function. At all orders, N remains 
a parameter. Because the N scaling part of the problem has been compartmentalized away 
from the rest of the physics, this work does not need to be repeated for different iV or for a 
different interaction potential. 

Also the derived Clebsch-Gordon coefficients allow this work for the correlated next- 
order wave function for systems with spherical symmetry to be extended tom = systems 
with cylindrical symmetry, such as most laboratory BECs (the additional normal modes 
associated with the z direction transform under the [N] and [N — 1, 1] irreps. of the Sn 
group). Furthermore, the fact that the Sn point group allows higher orders to be calculated 
also paves the way for an investigation of the above mentioned phenomenon of fermionization 
and crystallization [42], which display structure in the density profile which is absent from 
the correlated lowest-order density profile. Although, the present paper addresses L = 
confined quantum systems, in principle the formalism has been developed [43, 44] which 
would allow the extension of this work to higher angular momentum states. 

Systems of particles without a confining potential, but which bind under their own inter- 
action may be treated in a similar fashion, in Jacobi coordinates for example. 



VII. ACKNOWLEDGMENTS 

We gratefully acknowledge continued support from the Army Research Office. We thank 
David Kelle for advice on the use of graphs. 



37 



APPENDIX A: IDENTITIES 



+ IZu* ( Al ) 



dS 1 / 2 yv dv v d5 l l 2 



2 



d \ _, d d , d d ( d V 

d \ 3 d d d 

d^) -^^^5-5-5- ^ ^ (A3) 

_. , d d d _. d ( d \ ( d \ 

APPENDIX B: DERIVATION OF THE CLEBSCH-GORDON COEFFICIENTS 

Deriving the first Clebsch-Gordon coefficient C* 000 ' 7 is straightforward enough. It is simply 

C 000 ^ = 1 . (Bl) 

Strictly speaking, Clebsch-Gordon coefficients effect a unitary transformation and satisfy 
the normalization condition. 

We do not, however, need the above normalization relationship: we only need the coefficients 
to be linearly independent and to span the required space. Therefore, we drop the unitarity 
requirement and use unnormalized Clebsch-Gordon coefficients. 

Now according to the general theory of finite groups (Hamermesh, Section 5-4, p. 136), 
there is only one bilinear invariant of an irreducible representation. If this irrep. is unitary 
the invariant is proportional to the Kronecker delta function 5^ . This invariant couples 
two irreps. together to give a scalar [N] representation. 



Thus we have 



C^ 1 = S i3 (B3) 



38 



and 

C (ij){kl) = 6 ik$jl , (B4) 

where 1 < % < N - 1 and 1 < j < N - 1. 

More generally we can calculate the required Clebsch- Gordon coefficients by transform- 
ing the binary invariants B(Q) to symmetry coordinates with the required W matrices. We 
perform this and the remaining derivations in this Appendix symbolically using Mathe- 
matical]. 

Using the proportionality relation 

clS' 1 K WMwl] jm [w% n [B{g)] lmn , (B5) 

we define the unnormalized coefficient 

+t{t + i)® k -A,j + k(k + i)®j-A,k + ]{j + i)©*-Afc) , ( B6 ) 

where 1 < i < N - 1, 1 < j < N - 1, and 1 < k < N - 1. 
Similarly, using the proportionality relation 

CSS,* « [W 7 2 ] feV K 1 ] fen [W, 1 ]^[ J B(^)] mnp , (B7) 

we obtain the coefficient 



c: 



211,7 



^'^ ^ + l)(j-3)(j-2)A;(A; + l)/(/ + l) 
x (2i(i + l)(e_ j+ ,+i - e z _ fc )5 i)fe + 2i(i + l)(0_ j+fe+1 - e fc _,)5i,i 

+ i)(j - 2)(j - + - + i)e,_ fe 4,i 

+2*(* 2 - l)5 tAl ) , (B8) 

where 1 < i < j - 2, 4 < j < N, 1 < fc < N - 1, and 1 < I < N - 1 (therefore i < N — 2). 
The triple delta function tfj.j.fc = 1 when i = j = k and is zero otherwise. Note that the 
above expressions for the Clebsch-Gordon coefficients are valid for all N, not just a single 
N. In the context of this project, this is a very important point. 
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We could continue to derive the remaining coefficients in like manner, but due to the 
binary invariants and W transformations involved in these sectors, the increase in complexity 
is prohibitive. It is most convenient at this point to use the above Clebsch-Gordon coefficients 
to derive those that remain. We derive C? 2 l}{ n from 

(ij)(kl)m 

n 22\,I _ r <211,/ r <211,/ ^110,7^(111,7 /-r> Q \ 
^(ij)(kl)m ~ ^{ij)np^(kl)qr^pr ^nqm l Dy ) 

where all of the implied summations run from 1 to N — l and where 1 < % < j — 2, 4 < j < N , 
1 < fc < Z - 2, 4 < Z < A (therefore i<N-2k<N-2). 

The last two C 222 ' n Clebsch-Gordon coefficients C 222,1 , and C 222 ' 11 are a special case. 
Two linearly independent invariants may be formed as follows 

^222,7 _ ^211,7 ^211,7 ^211,7 ^110,7^110,7^110,7 / R1 n] 

(ii')(jj')(kk>) — U (ii')lm ^(jj')np U (kk')qr ^pq ^ rl \ Dl[J ) 

and 

^222,77 _ ^211,7 ^211,7 ^211,7 ^111,7^111,7 mm 

(H')(jj')(kk') ~ U (ii')lm U (jj')np ^ (kk')qr ^Inq ^mpr > \ Dll J 

where 1 < i < j - 2, 4 < j < N, 1 < k < I - 2, 4 < I < N,l < m < n - 2, 4 < n < N. The 
Clebsch-Gordon coefficients C 221 ' 1 , C 222,1 , C 222,n for arbitrary N have a more complicated 
expression than the others which we do not give here. A Mathematica[45] package for all 
the required Clebsch-Gordon coefficients may be found in Ref. 46. 



APPENDIX C: MULTIPLIERS (3(G) OF THE CLEBSCH-GORDAN COEFFI- 
CIENTS 

From Eq. (93), the transformed binary invariants Bw{Q) of Section VB are propor- 
tional to a Clebsch-Gordon coefficient of Sn in the irreducible symmetry coordinate basis 
of Section IV D (a linear sum of two Clebsch-Gordon coefficients in the 222 sector). In 
this Appendix, we derive these proportionality coefficients (3(G). These, together with the 
Clebsch-Gordon coefficients, will give us the symmetry coordinate transformed binary in- 
variants, B w (g) , of Eqs. (92) and (93). 
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By design, the transformations [Wj^,, have a simpler expression for smaller values of 
the row index v . Also, the multipliers (3(G) do not depend on the index v. Thus to 
derive /3 aia2Q ^(£) from Eq. (93), we use the lowest value for the index v in Eqs. (79), 
(80) (allowed by the row index restrictions on Eqs (62) , (64), (65)), along with the Clebsch- 
Gordon coefficients of Appendix B, which yield a non-trivial equation (two coupled equations 
in the case of the 222 sector - see below). 



1. Harmonic Order 



The only non-zero Clebsch-Gordon coefficients at harmonic-order which couple two or- 
thogonal irreps. of Sjy together to form a scalar [N] irrep. are proportional to 1, Ijv-i, and 
Ijv(jv-3)/2, and couple two [AT], [N — 1, 1] , and [N — 2, 2] respectively (I n is the nx n unit 
matrix). Thus we derive the proportionality coefficients from 



L J Al,A2 

Choosing the lowest value of v\, we obtain the results in Table IV. 



(CI) 



2. First Anharmonic Order 



The rank-one Clebsch-Gordon coefficient for a single [N] irrep. to give an [N] irrep. is 
obviously unity. Therefore, we obtain the proportionality coefficients from 

[W [ X N] UB(G)] V = \?P [N] (G)] „ , (C2) 



where Q is O or — • • From Eq. (C2) we obtain the results in Table V.] 

The rank-three Clebsch-Gordon coefficients have a significantly more complicated form, 



but the proportionality coefficients 



J Xl,X2,X3 

form using a computer algebra system such as Mathematica[45] and our own packages. 
We solve Eq. (93) for all but the 222 sector to give 



may still be computed in closed 



(1) ijaia2a3,I 



(G) 



Vl,V2,V3 



c; 



V\,V2,V 3 
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evaluated for the lowest z/ 1; u 2 , ^3 (allowed by the restrictions on the indices of W) for which 
the Clebsch-Gordan coefficient is non-zero. 



In the case of 
Eq. (93), i.e. 



(1) ^222,7 



and 



7,7,7 



(!) o222,77 



^222,77 ( £) 



, these may be evaluated from 



7.7.7 



(ij)(kl)(mn) 



(l)/o222,7 



^222,7 ^ [C 222 ' / ] fo -)(fc ( m „) + 



7,7,7 



(1) ^222,77 



<9) [c 



222,771 



7,7,7 



\(ij)(kl)(mn) 



(C3) 



for two different sets of ordered pair indices (ij)(kl)(mn) , where each of the ordered pairs 
are subject to the restrictions on the row indices of Eq (65). 



We list all the multipliers, 



, of the Clebsch-Gordan coefficients 



0) 00102013, TZ ^g'j 

J Xi,X2,X3 

in the expansion of the symmetry coordinate transformed, rank-three binary invariants, 



[Bw(G)\ 



V1,V2,V3 1 



in Tables VI-X. 



There is a subtlety regarding low values of N . When the number of vertices in a graph 
is less than N , the graph and corresponding binary invariant no longer exists (the vertices 
are associated with different values of the particle labels i j . . . , etc.) Since the transfor- 
mation from internal displacement coordinates to symmetry coordinates is a non-singular 
transformation, this reduction in the number of binary invariants is reflected in the number 
of non-zero, linearly independent Clebsch-Gordon coefficients at low N , i.e. the number is 
correspondingly reduced. These lower bounds on iV are noted in Tables VI-X. 

In the case of the expansion of the derivative portion of the kinetic term there is yet an- 
other subtlety; it is necessary to distinguish between graph edges corresponding to derivatives 
and those corresponding to non-derivatives. In this case, we use a vertical "tic" to mark 
derivative edges. There is only one graph which requires this distinction: 0~ . The binary 
invariant 5(0— ) is composed of the sum of binary invariants for the two graphs with two 
distinguishable edges: j3(0— ) = -8(0— ) + £?(0— ). Only the binary invariant B(<y*~) 
is present in the derivative term. The proportionality coefficient for this graph is shown in 
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TABLE I: Graphs for rank-two tensors at harmonic order, along with the corresponding tensor 
elements. 
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TABLE II: Graphs for rank-one tensors at first-order beyond harmonic, along with the correspond- 
ing tensor elements. 
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TABLE III: Graphs for rank-three tensors at first-order beyond harmonic, along with the corre- 
sponding tensor elements. 
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TABLE IV: Proportionality coefficients of the harmonic-order transformed binary invariants 
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TABLE V: Proportionality coefficients of the first anharmonic order transformed binary invariants 
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TABLE VI: Multipliers of the Clebsch-Gordon coefficients of the first-anharmonic-order trans- 
formed binary invariants: ^ paia 2 a 3 ,R ^g^ 
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TABLE VII: Multipliers of the Clebsch-Gordon coefficients of the first-anharmonic-order trans- 



formed binary invariants: 
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TABLE VIII: Proportionality coefficients of the first-anharmonic-order transformed binary invariants: 
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TABLE IX: Proportionality coefficients of the first-anharmonic-order transformed binary invari- 

, where ajPj = N\/(N - i)\ 
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TABLE X: Proportionality coefficients of the first-anharmonic-order transformed binary invariants 
with distinguishable edges: relevant graphs 
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